\chapter{AmberMDrun:自动化执行Amber MD的程序包}
\section{引言}
得益于高性能计算硬件的发展，分子动力学模拟越来越成为药物研发，微观机理探索等领域的重要方法。Amber \cite{amber}是一个程序集合的名称，它允许用户进行分子动力学模拟，特别是在生物分子方面。尽管没有单独的程序使用这个名字，但各个部分协同作用并提供强大的框架，使得许多常见的计算变得简便。此术语也应用于指代实现的经验力场。需要注意的是，代码和力场是分离的。除了使用Amber程序实现的力场外，还有其他计算机软件包也可以实施Amber力场。此外，力场是开源的，而代码则在授权协议下分发。

关于自动执行Amber程序，MOE提供了一个模块用于图形化设置Amber模拟和生成执行的shell脚本。Charmm-gui提供了更广泛适用的力场化网页界面，并提供简单的脚本用于运行Amber。此外，殷赋科技推出了网页版执行Amber MD的程序。然而，这些自动执行Amber MD的方法有些不够便利，扩展性也不强，并且某些方法过于昂贵。

因此，本章重点介绍AmberMDrun，这是一种高度可扩展和方便的程序套件，自动执行Amber MD模拟。它提供了方便易用的可执行程序，允许用户轻松计算蛋白质配体结合的自由能，并具有扩展性和便利性。
\section{理论方法}
Amber 分子动力学模拟软件包有两个主要的MD引擎。Sander是原始实现的版本，Pmemd是Sander的重构，主要提高代码性能，引入高性能GPU计算，实现性能飞跃。需要注意的是，除非有特殊的需求，否则大多数分子动力学模拟是在等温恒压（NPT）条件下进行的。
\subsection{Thermostat methods}
在MD模拟中，必须通过恒温器方法保持系统温度，其中通过调节粒子速度直接控制温度。在Berendsen恒温器中 \cite{doi:10.1063/1.448118}，温度由具有恒定温度的外部热浴控制，并且通过从热浴吸收或向热浴释放能量来维持系统的温度。当系统远离平衡时，温度调节非常有效，但粒子的动能并不严格遵循玻尔兹曼分布。因此，Berendsen恒温器可应用于加热阶段，但不适用于预平衡阶段。当系统完全平衡后，将使用弱耦合常数（例如10 ps），以便Berendsen恒温器可用于生产（采样）阶段。注意，Berendsen恒温器不适用于隐式溶剂模拟，因为它不能通过与溶剂分子的碰撞来帮助保持恒定的温度。

在Langevin恒温器中 \cite{doi:10.1063/1.1667473,doi:10.1021/ct800573m}，加热浴控制温度的速度和精度处于中等水平，它通过虚拟随机碰撞调整粒子的速度，扰乱系统的正常演化，削弱粒子之间的速度相关性。Langevin恒温器严格遵循规范集合属性，这不会影响每个状态的遍历性。因此，它适用于模拟热力学财产（如结合自由能），但不适用于模拟动力学财产。它在隐式溶剂计算中尤其有用，因为它可以模拟由溶剂引起的力。然而，当用显式溶剂运行全原子MD时，引入摩擦和随机力可能会对动力学产生不希望的影响。另一方面，由于作用的势能而获得“真实”动态的唯一绝对正确的方法是使用非常大的周期性水箱或更好的是非常大的水滴（以避免由于周期性而产生的伪影）在NVE模式下运行。虽然引入恒温器，如Langevin恒温器，可能会因势能而干扰运动，但通常有必要这样做，以在所需温度下准确模拟系统。根据所需的温度和模拟条件选择适当的摩擦力和随机力水平是至关重要的，这可能会因研究问题和所研究的系统而异。值得注意的是，正如Kleinerman等人所证明的那样，即使水摩擦减少了100倍，Langevin恒温器仍然可以产生准确的结果 \cite{10.1063/1.2943146}。

扩展拉格朗日方法可以利用时间反演克服上述方法的缺点。一种代表性的方法是Nosé-Hoover恒温器 \cite{PhysRevA.33.4253}，这是具有时间可逆性的最简单方法。然而，当系统远离平衡时，温度振荡很大，不容易收敛。因此，Nosé-Hhoover恒温器不适合加热阶段，但最适合平衡采样，因为它实现了规范系综的热力学，并且可以近似获得真实的动力学行为。值得一提的是，Nosé-Hoover恒温器在特定系统中可能有病态行为 \cite{doi:10.1021/ct8000365,doi:10.1021/ct400109a}。近年来，已经开发了一些改进的方法，例如优化的等速Nos-Hoover链（OIN）和随机等速Nos-胡佛RESPA积分方案 \cite{doi:10.1063/1.4848716,doi:10.1063/1.4999447}。
\subsection{Barostat methods}
在NPT系统中，有必要使用气压调节器来控制系统压力。主要使用的方法包括Berendsen、Nosé-Hoover和Parrinello-Rahman恒压器 \cite{doi:10.1063/1.3073889}。当系统远离平衡时，Berendsen恒压器在控制压力方面有效，适合系统的初始压力松弛，但不适合平衡采样，因为它无法实现标准分布的采样。与Berendsen方法相比，Nosé-Hhoover和Parrinello-Rahman方法更适合于平衡采样的压力控制 \cite{PhysRevA.33.4253,doi:10.1063/1.3073889}。蒙特卡洛（MC）方法也很有用，但其效率相对较低。因此，MC更适合于采样阶段 \cite{VITALIS200949}。
\subsection{十步平衡协议}
系统平衡由模拟准备的十个步骤组成 \cite{doi:10.1063/5.0013849}。MIN1是一个1000步最陡下降（SD）最小化，其中受体的重原子受到约束。MIN2是在15ps下的MD模拟，时间步长为1fs，对于大分子的重原子和涉及H原子的约束能为5.0kcal/mol。MIN3与MIN1相似，不同之处在于大分子的重原子结合力变为2.0kcal/mol。MIN4是一种1000步梯度下降能量优化，对大重原子的限制能量为0.1 kcal/mol。MIN5是在没有任何约束的情况下将1000个步骤的SD最小化。对于重原子和大分子中涉及H原子的键的能量为1.0 kcal/mol的NPT系统，MIN6是5ps MD，时间步长为1fs。MIN7类似于MIN6，不同之处在于大分子重原子上的能量为0.5千卡/摩尔。MIN8是MD模拟，NPT系综中的时间步长为1fs，时间为10ps。在MIN8中，蛋白质的非氢主链原子、核酸残基和大分子的重原子受到限制。MIN9是一个10ps的NPT模拟，时间步长为2fs，没有任何限制。MIN10是一种最终密度稳定的MD。
\subsection{MM-PB(GB)SA calculations}
分子力学/Pisson-Boltzmann（广义玻恩）表面积（MM/PB（GB）SA）方法通常用于计算非共价结合配合物的结合自由能 \cite{10.3389/fmolb.2017.00087}。配体和受体之间的结合自由能是配合物相对于受体和配体之间的自由能总和的相对自由能，通过公式 \ref{eq1}计算。Kollman等人开发了用于计算结合自由能的分子力学Poisson-Boltzmann表面积（MM-PBSA）和分子力学广义玻恩表面积的（MM-GBSA）方法 \cite{Srinivasan1998ContinuumSS,doi:10.1021/ar000033j,doi:10.1080/07391102.1998.10508279}。

MM-PBSA是计算不同分子体系自由能的有效方法 \cite{https://doi.org/10.1002/minf.201100135}，PBSA方法更耗时，而GBSA是PBSA方法的分析近似。溶剂化项的MM-PBSA和MM-GBSA自由能$\Delta G_{sol}^{complex}$, $\Delta G_{sol}^{receptor}$ ,$\Delta_{sol}^{ligand}$包括通过求解Poisson-Boltzman方程或广义Born方程计算的电子相互作用（极性贡献）和疏水相互作用（非极性贡献，\cref{eq1,eq2,eq3,eq4}）。
溶剂化相互作用包括电子（极性）相互作用$\Delta G_{GB/PB}$与静电势的相互作用，静电势通过求解Poisson Boltzman或广义Born方程计算，疏水（非极性）交互作用 $\Delta G_{non-polar}$ 通过溶剂可接触表面积（SASA）经验计算 \cite{https://doi.org/10.1107/S0021889883010985,Sitkoff1994AccurateCO,https://doi.org/10.1002/jcc.21372}。
	
其中$\gamma$和b是常数，可通过拟合实验烷烃转移自由能与SASA来确定。为了获得MM-PB（GB）SA结合自由能，可以使用单轨迹方案（STP）和多轨迹方案（MTP）两种方法，其中STP从复合物的轨迹中提取受体和配体轨迹，MTP使用复合物、受体和配体的多轨迹 \cite{doi:10.1021/acs.chemrev.9b00055,LEE2006864}。
	\label{eqall}
	\begin{equation}\label{eq1}
		\Delta G_{bind}=\Delta G_{bind,vacuum}+\Delta G_{sol}^{complex}- (\Delta G_{sol}^{receptor}+\Delta G_{sol}^{ligand}) 
	\end{equation}
	\begin{equation}\label{eq2}
		\Delta G_{sol}=\Delta G_{polar}+\Delta G_{non-popar}=\Delta G_{PB/GB}+\Delta_{non-popar}
	\end{equation}
	\begin{equation}\label{eq3}
		\Delta G_{non-popar}=\gamma \times \Delta SASA+b 
	\end{equation}
	\begin{equation}\label{eq4}
		\Delta G_{bind,vacuum} = \Delta E_{MM}- T\Delta S=\Delta E_{bonded} + \Delta E_{ele}+\Delta E_{vdW} - T\Delta S 
	\end{equation}
\subsection{AM1-BCC电荷计算方法}
AM1-BCC是一种准确而高效的方法，用于计算原子电荷，以表示分子的电子分布。AM1-BCC方法结合了Austin模型1(AM1)半经验方法生成初步猜测电荷，和键电荷增量 (BCC) 方案来修正这些电荷 \cite{model1,model2}。

在AM1-BCC方法中，特定原子（$j$）的原子电荷（$q_j$）是通过两个不同项的组合确定的（公式 \ref{eq5}）。 $q^{pre}_{j}$是通过快速预充电过程确定的初始值，该过程考虑了大多数化学方面，但仅使用该值不足以用于物质在其凝聚相中的模拟。AM1方法利用基于键对称性的平均原子电荷进行计算。然而，从AM1方法获得的静电电势（ESP）不够精确，无法与从更先进的方法（如\mbox{HF/6-31G$\ast$}获得的ESP相提并论。因此，在公式（\ref{eq5}）中添加了一个校正项$q^{corr}_{j}$，用于调整AM1原子电荷并接近重现HF/6-31G∗的静电势。 $q^{corr}_{j}$由公式 \ref{eq6}定义，其中$p\alpha$是键电荷校正（BCC），用于修复分子中各类型键的原子电荷。键连接模板矩阵T是用于计算从键电荷到以原子为中心的电荷的一种数学工具。它将分子中每种键的键电荷校正映射到相应的原子类型上。


总的来说，AM1-BCC是一种稳健的方法，可以为各种分子系统提供可靠的部分原子电荷。它被广泛用于药物发现和分子模拟，以检查分子与蛋白质和其他生物靶标的相互作用 \cite{doi:10.1063/5.0019056}。
\begin{equation}
	\label{eq5}
	q_j = q^{pre}_{j} + q^{corr}_{j} 
\end{equation}
\begin{equation}
	\label{eq6}	
	q^{corr}_{j} = \sum_{\alpha=1}^{\gamma} T_{j\alpha p\alpha}
\end{equation}
\section{结果与讨论}
我们程序的核心代码AmberMDrun是用C++编写的，通过Pybind11进行python绑定。用C++编写的核心代码确保了代码性能，python绑定为用户提供了易用性。代码为所有参数设置了合理的默认值。同时，用户可以通过修改一些参数来定制自己的MD模拟。对于希望扩展包的用户，良好的面向对象设计也使扩展成为可能。在下文中，我们将介绍我们程序的主要流程及其在MD模拟中的应用。
\vspace{-6pt}
\subsection{常用参数设置}
\begin{table}[htbp] 
	\caption{AmberMDrun中的一些关键参数。\label{tab1}}
	\newcolumntype{C}{>{\centering\arraybackslash}X}
	\begin{tabularx}{\textwidth}{CCCC}
		\toprule
		\multirow[m]{3}{*}{\textbf{Params}}	& \textbf{MIN}	& \textbf{NVT} & \textbf{NPT}\\
		& \textbf{imin = 1} & \textbf{imin = 0}&\textbf{ imin = 0} \\
		&&\textbf{ntb = 1} &\textbf{ntb = 2}  \\
		\midrule
		temp = 298.15 \textsuperscript{1}		& \ding{56}			& \ding{52} & \ding{52}\\
		cut = 8.0	\textsuperscript{2}	& \ding{52}			& \ding{52} & \ding{52}\\
		ntpr = 50	\textsuperscript{3}	& \ding{52}			& \ding{52} & \ding{52}\\
		ntwr = 500	\textsuperscript{4}	& \ding{52}			& \ding{52} & \ding{52}\\
		ntwx = 500 \textsuperscript{5}		& \ding{52}			& \ding{52} & \ding{52}\\
		maxcyc = 1000 \textsuperscript{6}		& \ding{52}			& \ding{56} & \ding{56}\\
		ncyc = 10	\textsuperscript{7}	& \ding{52}			& \ding{56} & \ding{56}\\
		ntmin = 10	\textsuperscript{8}	& \ding{52}			& \ding{56} & \ding{56}\\
		nstlim = 5000 \textsuperscript{9}		& \ding{56}			& \ding{52} & \ding{52}\\
		dt = 0.002	\textsuperscript{10}	& \ding{56}			& \ding{52} & \ding{52}\\
		irest = False \textsuperscript{11}		& \ding{56}			& \ding{52} & \ding{52}\\
		tautp = 1.0	\textsuperscript{12}	& \ding{56}			& \ding{52} & \ding{52}\\
		taup = 1.0	\textsuperscript{13}	& \ding{56}			& \ding{56} & \ding{52}\\
		gamma\_ln = 5.0	\textsuperscript{14}	& \ding{56}			& \ding{52} & \ding{52}\\
		nscm = 0 \textsuperscript{15}	 	& \ding{56}			& \ding{52} & \ding{52}\\
		ntc = 2	\textsuperscript{16}	& \ding{56}			& \ding{52} & \ding{52}\\
		ntf = 2	\textsuperscript{17}	& \ding{56}			& \ding{52} & \ding{52}\\
		thermostat \textsuperscript{18}		& \ding{56}			& \ding{52} & \ding{52}\\
		barostat \textsuperscript{19}		& \ding{56}			& \ding{56} & \ding{52}\\
		igamd = false \textsuperscript{20}		& \ding{56}			& \ding{56} & \ding{56}\\
		
		
		\bottomrule
	\end{tabularx}
	\noindent{\footnotesize{\textsuperscript{1} 模拟温度}}
	\noindent{\footnotesize{\textsuperscript{2} 指定非键截止值，单位为埃}}
	\noindent{\footnotesize{\textsuperscript{3} 以人类可读的形式打印到文件}}
	\noindent{\footnotesize{\textsuperscript{4} 重新启动文件输入频率。}}
	\noindent{\footnotesize{\textsuperscript{5} 坐标文件输入频率。}}
	\noindent{\footnotesize{\textsuperscript{6}最小化的最大循环次数。}}
	\noindent{\footnotesize{\textsuperscript{7} 在ncyc循环之后，最小化方法将从最陡下降转换为共轭梯度。 }}
	\noindent{\footnotesize{\textsuperscript{8} 最小化方法的标志。}}
	\noindent{\footnotesize{\textsuperscript{9} 要执行的MD步骤数。}}
	\noindent{\footnotesize{\textsuperscript{10} 时间步长。}}
	\noindent{\footnotesize{\textsuperscript{11} 用于重新启动模拟的标志。}}
	\noindent{\footnotesize{\textsuperscript{12} 时间常数。}}
	\noindent{\footnotesize{\textsuperscript{13}压力松弛时间。}}
	\noindent{\footnotesize{\textsuperscript{14} 碰撞频率。}}
	\noindent{\footnotesize{\textsuperscript{15} 以规则间隔移除平移和旋转质心（COM）运动的标志（默认值为1000）。}}
	\noindent{\footnotesize{\textsuperscript{16} SHAKE的标志，用于执行共价键的长度约束。}}
	\noindent{\footnotesize{\textsuperscript{17} 力评估。}}
	\noindent{\footnotesize{\textsuperscript{18} 恒温器选择。}}
	\noindent{\footnotesize{\textsuperscript{19} 恒牙器选择。}}
	\noindent{\footnotesize{\textsuperscript{20} 是否执行高斯加速动力学模拟。}}
\end{table}
\subsection{可扩展性}
开发人员很容易扩展AmberMDrun。大多数模拟是针对NPT系统的。我们使用高斯加速分子动力学（GaMD）增强采样来说明我们程序的易用性。我们从NPT类继承了GaMD类，以实现预期的GaMD模拟。首先，我们需要创建一个继承自NPT类的C++类GaMD。我们需要为GaMD唯一参数设置变量，例如igamd、ie等。然后我们需要做的是重写NPT中的公共和受保护的方法。对于大多数方法，我们只需要继承类NPT中的方法，而不进行更改。我们只需要在writeInput函数中写入GaMD添加的参数。需要做的另一件事是重写需要设置越来越多参数的构造函数。这就是为什么有必要通过python实现一个简单的调用，因为C++缺少关键字参数。最后，通过pybind11绑定到python，实现了类GaMD。可以看出，使用不同的采样方案执行MD模拟非常容易。
\subsection{AmberMDrun中的核心代码类}
\textbf{Systeminfo class}\quad 检查pmemd.cuda是否可用，然后使用pmemd.cud或sander生成运行Amber MD作业的命令。此外，在随后的系统平衡中，将获得受体残基的数量，用于受体的位置限制。

\textbf{MIN class} \quad 执行能量最小化。设置Maxcyc（最小化的最大循环数）和ntmin（最小化方法的标志）等（表 \ref{tab1}）.


\textbf{NVT class} \quad 执行NVT模拟。设置恒温器方法和nstlim（要执行的MD步骤数）等参数（表 \ref{tab1}）.


\textbf{NPT class} \quad 执行NPT模拟。设置气压计和恒温器方法以及nstlim（要执行的MD步骤数）等（表 
\ref{tab1}）.
\unskip
\subsection{基于AmberMDrun框架实现的计算MM-PB(GB)SA流程}
\begin{figure}[H]
	\centering
	\includegraphics[width=12cm]{figures/flow.pdf}
	\caption{AmberMDrun程序的工作流程，它依赖于AmberTools、Amber（可选）和gmx\_MMPBSA \cite{amber,doi:10.1021/acs.jctc.1c00645}。十步平衡协议被设置为默认的系统平衡方法，用户也可以定义自己的平衡方法。使用gmx\_MMPBSA、MM-PBSA或MM-GBSA结合自由能进行估算。\label{fig1}}
\end{figure}
我们通过编写简单的python代码，实现了一种自动化执行MM-PB(GB)SA计算的流程。通过使用AmberMDrun的库函数，整个流程代码仅仅需要几十行代码即可。首先，需要对蛋白和配体进行处理以符合Amber程序包的输入结构文件格式，因此，对蛋白调用pdb4amber，去掉氢原子，非标准残基以及水，重拍残基序号，这一步是必须的。对小分子则调用acpype进行力场分配，如果没有计算电荷，则调用sqm对小分子进行bcc电荷的计算。然后通过leap程序，将蛋白和配体组合在一起，输出Amber的拓扑和坐标文件。准备阶段完成之后，接下来就是需要进行平衡和生产项模拟，由于AmberMDrun提供的API，直接实现十步平衡法是非常容易的，仅仅需要实例化一些类即可，示例代码中提供了十步平衡协议的实现。因此，如果用户不想使用自己的平衡协议，通过默认的十步协议也是可以的。平衡过程结束之后，生产项模拟则仅仅需要实例化一个MD类，之后设定好参数就可以运行。生产模拟结束之后通过gmx\_MMPBSA软件包进行MM-PB(GB)SA结合自由能的计算。最后，对于模拟后的性质分析，也可以通过cpptraj的python接口pytraj进行，由于cpptraj的所有方法几乎全部绑定到python中，因此分析的代码同样是非常简单的。
\section{示例}
对于生物大分子研究和常用药物辅助设计，通常需要预测MM-PB(GB)SA。在这项工作中，AmberMDrun脚本工具提供了一个简单易行的脚本实现，从自动生成Amber输入、系统平衡到MD生产模拟，再到绑定自由能量计算。例如，Amber教程中的雷洛昔芬与雌激素受体结合的系统进行了测试。在AmberMDrun脚本中，只需要蛋白质和配体坐标文件。Amber输入文件是使用tleap和Acpype生成的。配体的力场参数化使用Acpype进行，配体使用gaff2的一般力场，当缺少小分子电荷时，通过bcc方法计算小分子电荷 \cite{https://doi.org/10.1002/jcc.20035,WANG2006247,model1,model2,doi:10.1063/5.0019056,SousadaSilva2012}。蛋白质的力场是amber14SB \cite{doi:10.1021/acs.jctc.5b00255}。系统平衡由（默认）十步协议自动执行 \cite{doi:10.1063/5.0013849}。然后，对生产进行50 ns MD模拟，最后，受体和配体之间的结合自由能是通过调用gmx\_MMPBSA获得的 \cite{doi:10.1021/acs.jctc.1c00645}。gmx\_MMPBSA最初旨在为Gromacs用户调用MMPBSA.py程序 \cite{doi:10.1021/ct300418h}提供方便。在我们的AmberMDrun脚本中，它被用于计算受体-配体复合物的结合自由能量。

\begin{figure}[H]
	
	\includegraphics[width=0.95\textwidth]{figures/MMPBSA_TOTAL.eps}
	
	\caption{使用AmberMDrun脚本，在raloxifen-雌激素受体的50纳秒分子动力学模拟中，比较使用Amber教程中的简单平衡方法所得到的MM-PBSA结果。\label{fig2}}  
\end{figure} 


% Example of a figure that spans the whole page width. The same concept works for tables, too.
\begin{figure}[htbp]
	
		\centering 
		
		%\includegraphics[trim=40 0 50 50,width=15cm]{Definitions/RMSD}
		%\includegraphics[trim=10 10 10 10,width=15cm]{Definitions/RMSF}
		\begin{minipage}[t]{0.5\linewidth}
			\centering		
			\includegraphics[width=\textwidth]{figures/FIG_RMSD}
			%\centerline{(a) RMSD}
		\end{minipage}%
		\begin{minipage}[t]{0.5\linewidth}
			\centering
			\includegraphics[width=\textwidth]{figures/FIG_RMSF}
			%\centerline{(b) RMSF}
		\end{minipage}
		
		

	\caption{Amber教程中的简单平衡方法和AmberMDrun脚本中的十步平衡协议的raloxifen-雌激素受体50 ns分子动力学模拟的RMSD和RMSF结果。\label{fig3}}
\end{figure}

在MD模拟中，通常使用RMSD（均方根偏差）和RMSF（均方根波动）来确定蛋白配体配合物的稳定性。根据图\ref{fig3}，使用十步平衡和简单平衡的雷洛昔芬雌激素受体复合物在整个模拟过程中是稳定的。使用十步平衡法的总体RMSD和RMSF通常低于Amber教程中雷洛昔芬雌激素受体示例的简单平衡 \cite{amber_tutorial}。


我们还在Amber教程中使用AmberMDrun脚本测试了雷洛昔芬雌激素受体复合物的例子。我们比较了表 \ref{tab2}中的三个结果:（1）来自本教程，（2）重复Amber教程中的模拟平衡准备，以及（3）十步平衡\cite{doi:10.1063/5.0013849}。可以看出，结果（2）和（3）几乎相同，结果（1）的结合自由能略高于（2）或（3）。使用脚本工具AmberMDrun，只需单击一次即可完成复杂的动力学模拟 \verb|mmpbsa -p protein.pdb -m mol.mol2 --ns 50|，它是自动执行的，实现了我们轻松运行Amber的目标，就像openmm一样，它可以通过Python实现其功能 \cite{doi:10.1021/ct300857j}。
\begin{table}[H] 
	\caption{MM-PBSA和来自雷洛昔芬-雌激素受体复合物的2ns MD模拟的组分结果的比较。\label{tab2}}
		\newcolumntype{C}{>{\centering\arraybackslash}X}
		\begin{tabularx}{{\textwidth}}{clllllll}
			\toprule
			\textbf{Result}	& \boldmath\textbf{$\Delta vdW$}	& \boldmath${\Delta E_{EL}}$ & \boldmath${\Delta E_{PB}}$ & \boldmath${\Delta E_{NPOLAR}}$ & \boldmath${\Delta E_{GAS}}$ & \boldmath${\Delta E_{SOL}}$ &  \textbf{$\Delta$TOTAL} \\
			
			\midrule
			教程 \textsuperscript{1}&$-$56.29&$-$38.03&57.64&$-$5.04&$-$94.32&52.60&$-$41.73\\
			简单平衡 \textsuperscript{2} &$-$55.48&$-$41.00&53.79&$-$5.58&$-$96.48&48.21&$-$48.27\\
			十步平衡协议 \textsuperscript{3} & $-$58.62&$-$36.00&54.00&$-$5.70&$-$94.63&48.30&$-$46.32\\
			\bottomrule 
		\end{tabularx}

	\noindent{\footnotesize{\textsuperscript{1} Amber官网教程中的结果。}}
	\noindent{\footnotesize{\textsuperscript{2} 使用教程中平衡输入文件的结果。}}
	\noindent{\footnotesize{\textsuperscript{3} 使用十步平衡协议的结果。}} 
	
	
\end{table}
\section{扩展代码实现自己的类}
我们将结合伪代码和C++代码，通过高斯加速分子动力学说明如何扩展AmberMDrun自定义类，实现特殊的功能。首先，高斯加速分子动力学在Amber中的关键字为igamd，表示是否启用高斯加速分子动力学模拟。因此，我们需要添加该参数，首先继承NPT类，并且实现构造函数的时候，采用C++中的默认参数。然后，为各个参数设置Get Set函数，所有的Get Set方法均可以返回动态的类型转换，直接调用父类的方法，除非有特殊的需求。并且重写父类中的wirteInput方法，剩余的函数均可以使用父类的方法。详细的代码在附录中，我们提供了gamd.hpp和gamd.cpp，可以看到构造函数的参数非常多，对于C++来说，并没有类似Python的关键字参数，所以使用C++的构造函数需要注意参数的顺序，当然也可以使用少量的初始化值，然后使用Set方法设置自己想要的参数，但是显然不如Pytho里的关键字参数便捷。因此，我们提供了绑定Pyhton的接口，借助于C++模板编程的强大的表达能力，为C++绑定到Python也变得非常简单。这三个类实现之后，用户自己扩展的类变完成了，对于用户来说几乎不需要编写很复杂的代码，也不需要很多计算机的知识。

总的来说，实现自己的类只需要三步：

1.实现自定义类的头文件，定义子类和子类中的方法。

2.为该自定义类实现源文件，主要是各个方法的具体实现。

3.通过Pybind11将其绑定到Python。

\section{本章小结}
AmberMDrun是一个由Python和C++编写的集成脚本，旨在通过处理蛋白质和小分子的结构、产生小分子力场、平衡系统、运行MD模拟和实现MM-PB（GB）SA结合自由能来自动生成输入文件。这个脚本是一个轻便方便的包，可以很容易地融入分子模拟研究，以加快进展。


对于增强型采样，例如GaMD，用户可以从现有的NPT类继承并实现自己的GaMD类。然而，许多用户不熟悉Python或C++，这可能是使用的障碍。未来，一个研究方向可能是为各种增强的采样动力学开发独特的类别。另一种选择可以是简化参数或引入图形界面，使用户更容易访问。此外，用C++重新编写tleap的主要组件，可以利用Python提供的用户友好的语法和图形界面，为生成拓扑坐标文件提供更方便的工具。此外，使用web服务器而不是图形界面也是一个不错的选择，正如文章 \cite{article}展现的形式。
